#include  <stdio.h>
#include  <math.h>
#include  <stdlib.h>
#include <sys/io.h>
#include <unistd.h>
#include <string.h>
#include <stdbool.h>
//#include <syslib.h>
//#define Npix  3.2   // 1 arcsec
//#define Notpix 44 //120  //60, 100*0.52 arcsec. For mini-GWAC, it corresponds to one pixel  error
#define Notpix 120 //120  //60, 100*0.52 arcsec. For mini-GWAC, it corresponds to one pixel  error
main()
{
    int     i,j,k,m,N,N1,N2,N3,F;
	float	rxc1,ryc2,rRA3,rDEC4,rB15,rB26,rR17,rR28,rI9;
	float   xc1,yc2,mag3,merr4,elli5,class6,fwhm7;
    float   otxc1,otyc2,otRA3,otDEC4;
    char    ottype5[200];
    int     flag8;
	float   delx,dely,delr,delotx,deloty,delotr;
    float   Npix;
    float   sum_B2,sum_R2,sum_comp,sum_I;
    float   average_comp,average_R2,average_B2,average_I;
    float   x[50000], y[50000],mag[50000],merr[50000],elli[50000],class[50000],fwhm[50000];
    int     flag[50000];
    float   ref1[50000],ref2[50000],ref3[50000],ref4[50000],ref5[50000],ref6[50000],ref7[50000],ref8[50000],ref9[50000];
    float   otx[1],oty[1],otRA[1],otDEC[1];
    char    ottype[100][200];
    FILE    *fp1,*fmatch,*fnomatch,*fcomp,*fdeltamag;
	char    *refall,*newimg,*otcoord;

	refall="refusno.cat";
	newimg="newimage.cat";
    otcoord="otcoord.cat";

//	printf("refall = %s\n",refall);
//	printf("newimg = %s\n",newimg);
//  printf("otcoord = %s\n",otcoord);

	i=j=k=m=N=N1=N2=N3=F=0;
    sum_comp=sum_R2=sum_B2=sum_I=0.0;
//    sum_delmag_R2=sum_delmag_B2=0.0;	
	fmatch=fopen("match_c.db","w+");
    fcomp=fopen("comp_c.db","w+");
	fnomatch=fopen("nomatch_c.db","w+");
    fdeltamag=fopen("deltamag_c.db","w+");
	fp1=fopen(newimg,"r");	
	if(fp1)
	{
	        while((fscanf(fp1,"%f %f %f %f %f %f %f %d\n",&xc1,&yc2,&mag3,&merr4,&elli5, &class6, &fwhm7, &flag8))!=EOF)
	        {
                x[i]=xc1;
                y[i]=yc2;
                mag[i]=mag3;
                merr[i]=merr4;
                elli[i]=elli5;
                class[i]=class6;
                fwhm[i]=fwhm7;
                flag[i]=flag8;
		        i++;
		    }
		N1=i;
	}
	fclose(fp1);
	
	fp1=fopen(otcoord,"r");	
	if(fp1)
	{
	        while((fscanf(fp1,"%f %f %f %f %s\n",&otxc1,&otyc2,&otRA3,&otDEC4,ottype5))!=EOF)
	        {
                otx[1]=otxc1;
                oty[1]=otyc2;
                otRA[1]=otRA3;
                otDEC[1]=otDEC4;
		        sprintf(ottype[1],"%s",ottype5);
            }
	}
	fclose(fp1);
	
    //	 refall is for the usno catalog
	fp1=fopen(refall,"r");
        if(fp1)
        {
          while((fscanf(fp1,"%f %f %f %f %f %f %f %f %f\n",&rxc1,&ryc2,&rRA3,&rDEC4,&rB15,&rB26,&rR17,&rR28,&rI9))!=EOF)
          {
          ref1[j]=rxc1;
          ref2[j]=ryc2;
          ref3[j]=rRA3;
	  ref4[j]=rDEC4;	
	  ref5[j]=rB15;
	  ref6[j]=rB26;
	  ref7[j]=rR17;
	  ref8[j]=rR28;
          ref9[j]=rI9;
          j++;
          }
          N2=j;
        }
        fclose(fp1);

	for(i=0;i<N1;i++)
	{
        delotx=x[i]-otx[1];
        deloty=y[i]-oty[1];
		delotr=sqrt(pow(delotx,2)+pow(deloty,2));
		k=0;
	    for(j=0;j<N2;j++)
		{
			delx=x[i]-ref1[j];
			dely=y[i]-ref2[j];
			delr=sqrt(pow(delx,2)+pow(dely,2));
			if(ref8[j]<8.0)
			{
				Npix=7.0;
			}
			else
			{
				Npix=3.2;		
			}
        //    printf("===========x=%.3f y=%.3f  delr=%.3f  i=%d\n",x[i],y[i],delr,i);
            		if(delr<Npix && delotr>Notpix)  //the obj in new image matched an obj in the usno but not in the error circular of miniOT
			{
				fprintf(fmatch,"%.3f %.3f %.3f %.3f %.3f %.3f %.3f %d %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f\n",x[i],y[i],mag[i],merr[i],elli[i],class[i],fwhm[i],flag[i],ref1[j],ref2[j],ref3[j],ref4[j],ref5[j],ref6[j],ref7[j],ref8[j],ref9[j],delr);
				k++;
                if(merr[i]<0.1 && merr[i]>0.001 && x[i]>50 && x[i]<1950 && y[i]>50 && y[i]<1950)  //to get the comp stars, which should be not fait objects and near the edges of the image.
                {
                    N++;
                    sum_comp=sum_comp+mag[i];
                    sum_R2=sum_R2+ref8[j];
                    sum_B2=sum_B2+ref6[j];
                    sum_I=sum_I+ref9[j];
                    fprintf(fcomp,"%.3f %.3f %.3f %.3f %.3f %.3f %.3f %d %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f\n",x[i],y[i],mag[i],merr[i],elli[i],class[i],     fwhm[i],flag[i],ref1[j],ref2[j],ref3[j],ref4[j],ref5[j],ref6[j],ref7[j],ref8[j],ref9[j],delr);
		//		    printf("=====%.3f %.3f %.3f %.3f %.3f %.3f %.3f %d %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f\n",x[i],y[i],mag[i],merr[i],elli[i],class[i],fwhm[i],flag[i],ref1[j],ref2[j],ref3[j],ref4[j],ref5[j],ref6[j],ref7[j],ref8[j],ref9[j],delr);
                }
			}
            if(delr<Npix && delotr<Notpix) //this obj is in the usno and in the error circular of miniOT
            {
                F++;
                k++;
                //printf("CATAS  Npix=%.3f  Notpix=%.3f\n",Npix,Notpix);
			    fprintf(fnomatch,"%.3f %.3f %.3f %.3f %.3f %.3f %.3f %d %.3f %.3f %.3f %s %d\n",x[i],y[i],mag[i],merr[i],elli[i],class[i],fwhm[i],flag[i],ref6[j],ref8[j],ref9[j],"CATAS",F);	
            //CATAS means that the obj is at the known catalog 
            }
		}
		if(k==0 && delotr<Notpix)  //this obj is not in the usno and in the error circular of miniOT
	    {            
            F++;
                //printf("Mini  Npix=%.3f  Notpix=%.3f\n",Npix,Notpix);
	    	fprintf(fnomatch,"%.3f %.3f %.3f %.3f %.3f %.3f %.3f %d %.3f %.3f %.3f %s %d\n",x[i],y[i],mag[i],merr[i],elli[i],class[i],fwhm[i],flag[i],99.000,99.000,99.000,"miniOT",F);	
        }
        if(k==0 && delotr>Notpix)  //this obj is not in the usno and also not in the error circular of miniOT,might be a new OT or a moving obj.
		{
            F++;
            fprintf(fnomatch,"%.3f %.3f %.3f %.3f %.3f %.3f %.3f %d %.3f %.3f %.3f %s %d\n",x[i],y[i],mag[i],merr[i],elli[i],class[i],fwhm[i],flag[i],99.000,99.000,99.000,"newOT",F);	
	    }
    }

    if(N>0)  //have comp stars
    {
        average_comp=sum_comp/N;
        average_R2=sum_R2/N;
        average_B2=sum_B2/N;
        average_I=sum_I/N;
	    fprintf(fdeltamag,"%.3f  %.3f %.3f %.3f %d\n",average_comp, average_B2, average_R2, average_I, N);	
    }
    else  //no comp stars
    {
        average_comp=average_B2=average_R2=average_I=0.0;
	    fprintf(fdeltamag,"%.3f  %.3f %.3f %.3f %d\n",average_comp, average_B2, average_R2, average_I, N);	
    }
    	    
	fclose(fmatch);
	fclose(fnomatch);
    fclose(fcomp);
    fclose(fdeltamag);


}
